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ABSTRACT 

We consider a two-parametric family of radially anisotropic models with non¬ 
singular density distribution in the centre. If highly eccentric orbits are locked near 
the centre, the characteristic growth rate of the instability is much less than the 
Jeans and dynamic frequencies of the stars (slow modes). The instability occurs only 
for even spherical harmonics and the perturbations are purely growing (aperiodic). 
On the contrary, if all orbits nearly reach the outer radius of the sphere, both even 
and odd harmonics are unstable. Unstable odd modes oscillate having characteristic 
frequencies of the order of the dynamical frequencies (fast modes). Unstable even 
harmonics contain a single aperiodic mode and several oscillatory modes, the aperiodic 
mode being the most unstable. 

The question of the nature of the radial orbit instability (ROI) is revisited. Two 
main interpretations of ROI were suggested in the literature. The first one refers to 
the classical Jeans instability associated with the lack of velocity dispersion of stars 
in the transverse direction. The second one refers to Lynden-Bell’s orbital approach 
to bar formation in disc galaxies, which implies slowness and bi-symmetry of the 
perturbation. Oscillatory modes, odd spherical harmonics modes, and non-slow modes 
found in one of the models show that the orbital interpretation is not the only possible. 

Key words: Galaxy: centre, galaxies: kinematics and dynamics. 


1 INTRODUCTION 

Spherical systems with predominance of eccentric orbits are 
subject to the so-called radial orbit instability (ROI) that 
leads to formation of non-spherical structures. These struc¬ 
tures are naturally associated with the triaxial bulges and 
bars, observed in a variety of self-gravitating systems, where 
the eccentric orbits could occur as a result of the radial col¬ 
lapse in early stages of formation. Numerical simulations 
of collapsing systems have been carried out in Aguilar and 
Merritt (1990), Roy and Perez (2004), Trenti and Bertin 
(2006). In addition to these non-equilibrium systems, it is 
of interest to study the stability of equilibrium models that 
are used in modelling galaxies, globular and open clusters. 
The stability conditions may impose substantial restrictions 
on the allowed parameters of models. 

By analogy with the well-known Ostriker - Peebles sta¬ 
bility criterion for disk systems, Polyachenko and Shukhman 
(1981) proposed a global stability criterion for spherical sys- 
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terns, Gut = 2 T r /T± > 1.7 ± 0.25, where T r and T± are 
total kinetic energy of radial and transverse motion. Sub¬ 
sequently it was found that the specified range for Grit is 
not rigorously found. In particular, a critical value of £ for 
generalized polytropic models was found close to 1.4 (Frid¬ 
man and Polyachenko, 1984, Barnes et al. 1986), or even 1 
(Palmer and Papaloizou, 1987)0 On the other hand, Osip- 
kov Merritt models give examples of systems that preserve 
initial spherical shape with Grit as much as 2.5 (Meza and 
Zamorano 1997). The most stable radially anisotropic con¬ 
figuration (Grit ss 2.9) were obtained by Trenti and Bertin 
(2006) in the numerical simulation of collisionless collapse. 
In the latest models, the stabilizing effect was due to the 
nearly isotropic core, while large anisotropy was achieved 
due to a strongly anisotropic shell. 

Several mechanisms were proposed to explain ROI, 
among which we mention two. The first mechanism treats 
this instability as the Jeans instability of anisotropic 


1 Note that the result by Palmer and Papaloizou was questioned 
by Polyachenko et al. (2011). 
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medium, in which the velocity dispersion in the transverse 
direction cannot resist gravitational attraction (Polyachenko 
and Shukhman, 1972, 1981; Antonov, 1973, Barnes et al. 
1986). The other mechanism supported by Merritt (1987, 
1999), Saha (1991), Weinberg (1991), Palmer (1994), and 
others is claimed to be similar to bar formation in rotating 
discs described by Lynden-Bell (1979): this is a tendency of 
some orbits to line up in the direction of the bar that leads 
to increase of the density and potential perturbation. 

A detailed exposition of the latter approach is given by 
Palmer (1994), who describes allowed mathematical simpli¬ 
fications of a general matrix equation for eigenoscillations to 
study instability under question. It was assumed that: 

• The modes are even (even spherical harmonics Z), since 
the stellar orbits are symmetric with respect to the centre 
of the system. 

• The modes are ‘slow’, i.e. the modes eigenfrequency u> 
must be much smaller than the characteristic radial fre¬ 
quency of stars, |a>| <C fli. Under this condition, periods 
of stars will be much shorter than characteristic time of the 
instability, and the orbits can be regarded as separate ob¬ 
jects. 

These two features certainly narrow the range of possible 
unstable modes compared to the more general former ap¬ 
proach, in which the frequencies w can be of the order of the 
radial frequency, |<n| ~ fli. For convenience, we denote the 
two mechanisms as ‘Jeans’ and ‘Lynden-Bell’, but the pri¬ 
mary difference is fastness and slowness of unstable modes. 
Besides the mentioned above restrictions, the ‘Lynden-Bell’ 
mechanism assumes that a radial part x( r ) of the eigenfunc¬ 
tion of the perturbed potential <5<f>(r, 9, tp) = x( r ) Y™(9, tp) 
is nodeless. 

In what follows we analyse the spectra of radially- 
anisotropic DFs of the form F(E, L) to answer the question 
whether only slow modes, or both types of modes are pos¬ 
sible. We imply that all unstable non-radial modes are due 
to ROI, provided they are stabilized by decreasing the ra¬ 
dial anisotropy. For the analysis we employ a two-parametric 
family of radially anisotropic models without central singu¬ 
larity (Polyachenko et al. 2013): 

F (E, L) = N ^ 3 L J t ] H ( L t ~ L)Fo(E) , 

F 0 (E) = 2(1+ q)(-2E) q , (1.1) 

where H(x) is the Heaviside step function; parameters q > 
— 1, Lt ^ 0; E is the energy and L is the absolute value of 
the angular momentum of individual stars, 

Em i(tv +v 2 x) + $ 0 (r) , L = rv±. 

In the distribution function (DF) (11.111 . the gravitational 
potential 'FoM is set to zero on the sphere boundary; 
the boundary radius R, and the total mass M are unity; 
N(q, Lt) is the normalization constant. 

In contrast to the well-known Osipkov - Merritt models 
(Osipkov, 1979; Merritt, 1985), and generalized polytropic 
models (Cumin, 1952) 

F GP (E,L)=C(s,q)L- s (-2E) q , (1.2) 

where —1 ^ q ^ 7/2, — oo < s < 2, our family of models 
lull have two advantages needed for correct analysis of the 


ROI: (i) the central density and the potential are finite, and 
(ii) for a wide range of parameter q, variation of parameter 
Lt transforms the system from isotropic to purely radial. 
A technical advantage is that the specific form of the In¬ 
dependence of the DF significantly simplifies a cumbersome 
numerical procedure of finding the eigenmodes. 

Below we investigate the stability of two series of DF 
HU with fixed values of the parameter q = 0 and q = —1. A 
principal difference between the two is in the energy distri¬ 
bution of stars: the former one holds equipartition, while 
the latter is mono-energetic, since lim Fo(E) = 5(E) 

9 - 1 - 1 + 

(see, e.g., Gelfand and Shilov, 1964). For highly radially 
anisotropic systems, models of q = 0 series have orbits with 
small apocentric distances (short needles) confined in the 
centre, where the characteristic orbital frequency fldyn is 
very large. In contrast, in mono-energetic models length of 
all highly eccentric orbits is nearly equal, and all stars can 
reach the outer radius of the system. We find that this fea¬ 
ture results in a completely different character of the insta¬ 
bility. 

The plan of the paper is as follows. In Section 2 we 
provide a general matrix equation to determine the eigen- 
frequencies, and its special form for the series under con¬ 
sideration. In Section 3 we present results of calculations of 
eigenmodes. Finally, in section 4 we summarize the results 
and discuss the physical mechanisms of radial orbit instabil¬ 
ity. 


2 MATRIX EQUATIONS 

In this paper, we address the problem of stability by finding 
the eigenfunctions x( r ) and the eigenfrequencies to of colli¬ 
sionless systems using matrix method. The method was first 
proposed by Kalnajs (1977) for disk systems. For spherical 
systems we are interested in a similar matrix equation, which 
was first obtained by Polyachenko and Shukhman (1981). 
Details of the derivation has been repeatedly given in the lit¬ 
erature (see, e.g., Polyachenko and Shukhman, 1981; Wein¬ 
berg, 1991; Bertin et al., 1994; Saha, 1991; Palmer, 1994), so 
here we only present the equation and explain the notations. 
The equation can be written as 

Det ||5 a/3 — A4“^(u;)|| = 0 , 3 =1,2,3,..., (2.1) 

where matrix A4 a ^(co) is given by the following expression 


OO Z /> 

M ap (u) = -4nG(2nf E E ^ / / 

Zi =—oo Zo =—l 


dEdL 


nE,L)[n hh A + ,^\( 1 ^-). (2.2) 


d 




id — f2l 1 i 2 


The integration in (12.211 is taken over the allowed domain 
T> of two-dimensional action sub-space (E,L). This is the 
so-called ‘Lagrangian’ form of the matrix elements (see Ap¬ 
pendix A). It is different from the more familiar ‘Euler’ form: 


M af3 (uj) = 4n G (2n) z 


oo Z /i 

E E A" / / 

— - 7 _ 7 J J 


Zi = — oo Zo = — Z 


dELdL 
fii (E,L)‘ 


T .dF , dF 

^Zl l 2 (E,L) + I2 

K h (E, L) - df dL , (2.3) 

« - ^hi 2 (E,L) 
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which can be formally obtained by integrating by parts and 
discarding the boundary terms. Note that ‘Euler’ form be¬ 
comes incorrect for F(E, L) with an integrable singularity 
at the sub-space boundary. For our series of models, this 
occurs for q < 0 at the boundary E = 0. 

In equations (12.211 and (12.311 indeces a and /3 correspond 
to the expansion of the radial part of the potential x( r ) 
and the radial part of the perturbed density II(r) over the 
biorthogonal set: 

x(r) = ^2C a x a {r) , (2.4) 

a 

n ( r ) = ^E c, VM ' (2 ' 5) 


The perturbations are assumed to be independent of 
the azimuthal variable ip, since the eigenfrequencies of 
the perturbations u> for spherically symmetric distribu¬ 
tions F = F(E,L) are independent of the azimuthal 
number m. Therefore, instead of the angular dependence 
of the general form <54? = x( r ) Fi m (0, +) e~ iut , 5p = 
II(r) Yj m (0, ip) e~ lu,t , one can consider simplified axisymmet- 
ric one <5$ = X (r) Pi(cosd) Sp = II(r) P;(cos6>) e~ iut . 

Here, functions ,\ a ( r ) and p a (r) are related by the Poisson 
equation 


d 2 2 d_ 
dr 2 r dr 


1(1 + 1) 


and satisfy the so-called biorthonormal conditions, 
i 

(r) p 13 (r) r 2 dr = —<5“^ . 

Note that x* and p a depend on the index l, but we omit it 
for brevity. 

Subscripts h and Z 2 correspond to the decomposition 


<54>(7i,7 2 , 101 , 11 ) 2 ) = 2 (I) exp[i(Ziwi + I 2 W 2 )] 

I 1 I 2 

of the spatial dependence of the perturbed potential in har¬ 
monics of angular variables wi and u >2 conjugate to action 
variables 7i and 7 2 : 

r max j 

h = ^f PrdT = i / \J 2E - 2 Mr)-^dr, I 2 = L . 

r min 


The radial angular variable wi is related to the radius r as 
follows: 

wi = fii 

An explicit expression for the angular variable W 2 can be 
found in the mentioned above papers (e.g., Polyachenko & 
Shukhman, 1981). For the perturbations independent of the 
azimuthal variable <p , <54? and <5 p do not depend on the an¬ 
gular variable W 3 . Functions Lli x i 2 (E,L) denote linear com¬ 
binations of orbital frequencies, = iifli + I 2 LI 2 , which 

are determined by: 

T max 

117 dr 

ni n r l j2E-2$ 0 (r)-L 2 /r 2 ’ 

' min V 


dr' 


y / 2E-24>o(r)-L 2 /r' 2 


and, 


' ma: 

02 _ l r 

fll TV J 


dr 


_ A p 


r 2 ,/2£-24>o(r)-75 2 /r 2 


where Aip(E, L) is the angular distance in the orbital plane 
between r m i n and r max . For highly eccentric orbits and non¬ 
singular unperturbed potentials the cases in which we are 
interested in, this angle is close to 7 7r, giving the frequency 


ratio n 2 /fii 


;. Such orbits, called 2:l-orbits, are slowly 


precessing ellipses symmetric relative to the centre. 

The coefficients Df are nonzero only for even |Z — fc| and 
equal 

D k = j_ (i + mi - fc)! 

22 ' [da-fc))!(la + ^)) ! ] 2 

Finally, iP?f h (E,L) = <f>? lh (E, L)<tf lh (E, L), where, 

7T 

<A“i i 2 ( E > L ) =* ^ J cosOi 1 i 2 (E,L-,wi)x°‘[r(E,L,wi)] dwi , 

0 

and the angle 0z 1 z 2 (E, L , wi) is defined by 0^ i 2 (E, L\ wi) = 

W\ 

fiqz 2 - hSip(E, L\ Wi), where 


r (E,L,wi ) 


5ip(E,L,wi) = L J - j 

7-mi n(E,L) X \ 


dx 


[2E + 24/(a;)] x 2 — L 2 


is the angular distance between r m i n and the current radius 
r; the relative potential 4/(r) = —4>o(r) > 0. 


3 RESULTS 

Equilibrium models of q = 0 and q = — 1 were analysed in 
detail in Polyachenko et al. (2013). In both cases, limiting 
models Lt = 0 describe systems of purely radial orbits with 
global anisotropy £ = 1 — (j -1 = 1. Nearly radial models 
corresponding to small Lt have isotropic and almost homo¬ 
geneous kernel with radius r\ ~ Lt- Within this radius, the 
potential is almost constant: 'l'(r) = 4/(0 ) + 0(r 2 ), 4/(0) > 0. 
For r > n, the density and the potential vary as follows: 

4/9+ 1 / 2 

^— , 4> ~ ln"(l/r) , n = (1/2 — q) . 

r 2 


3.1 q = 0 series 

The DFs F(E,L) in q = 0 series are independent of the 
energy E and the angular momentum L within the allowed 
domain X>: 


F(E,L) 


77(0, L t ) H(L t - L) 

2^ L 2 . 


H(—2E). 


(3.1) 


For the chosen form of DF, the expression for ma¬ 
trix elements M“ /3 (oj) is particularly simple, since the 
two-dimensional integration over T> is reduced to one¬ 
dimensional integration along two boundary lines: the ver¬ 
tical 0 < L < Lt, E = 0, and the horizontal L = Lt, 
E c < E < 0 (shown by thick lines in Fig. |TJ. The models 
become isotropic when the parameter Lt ^ 0.6682 (Poly¬ 
achenko et al., 2013). 
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Figure 1. ( E,L ) sub-space and the allowed domain T> (filled 
area). The upper curve is the line of circular orbits. 


A suitable expression for the matrix elements M aP (uj) 
can be obtained from (12.31) and written in a form containing 
eigenfrequency squared and summation over non-negative l\ 
only: 


M ap (uj) 


[ LdLS aP (0,L ) 
Lrp J 
0 


167V 

Lt 




JL L fr ,;r ^<^(0^) 

L T y f2r(0,L) [w 2 -0 2 ; (0,L)] 

0 1 2 


0 



E c 


n hl2 (E,L T )tff l2 (E,L T ) \ 
Q 1 (E,L t ) [w 2 - Qf ih (E,L T )] J 


(3.2) 


Here, the coefficients S'q =1/2 for Zi = 0, and Sq = 1 
otherwise. The first term on the r.h.s. of (IQIl results from 
a summation of terms independent of uo 2 : 


S ap (E,L) 


2 

fit (E,L) 


OO l 


J2s h D\^r p 2 {E,L), 


which can be performed analytically (see, e.g., Saha 1991): 


i r r 

* P (E,L) = - [ 


X“(r)x /3 (r) 


[2E — 2<E> 0 (r) — L 2 /r 2 ] 1 / 2 


dr . (3.3) 


Using equations m with matrix elements (El , we 
investigate the stability of spherical harmonics l in the range 
1 ^ l ^ 20 for Lt -S 0.01. The model Lt = 0.01 is highly 
radially anisotropic with ratio / = 2 T r /T± = 42, and the 
global anisotropy £ = 0.98. 

Success of the matrix method depends largely on the ap¬ 
propriate choice of basis functions {x“, p a }. In Appendix B 
we describe a method for constructing a variety of basis sets 
in which p a = A a g(r)x a , where function g(r) is arbitrary. 
The special case g(r) = — 1 corresponds to a well-known or¬ 
thogonal system of spherical Bessel functions (Polyachenko 
and Shukhman, 1981). 

A numerical code for mode’s calculation was tested by 
finding lopsided l = 1 shear modes, which is present in all 
models. They are called zero modes, since its eigenfrequency 
uj = 0, and the eigenfunction x = X' 1 ( r ) = Avk'. Recall 


q=0 series 



Figure 2. The dependence of the growth rate 7 of aperiodic 
modes on the parameter Lt for l = 2, 4, 6 . Stabilization of the 
most unstable mode 1 = 2 occurs at L cr ; t = 0.316. Profiles of 
characteristic frequencies are given by ( 13.51) . The shaded area 
shows the range of variation of the radial frequency Hi, (ffy) 1 . 0 c■ 
which almost coincides with the maximum radial frequency Hi 
(not shown). Results of the growth rate TV-body calculation for 
l = 2 are given by asterisks. 


that shift 5z of the sphere as a whole along z-axis gener¬ 
ates a perturbation of the potential <5$ = 5z ■ cos(9<l>o = 
5z <I>o Pi (cos 6). For small a; 2 , the matrix elements can be 
represented as a series in uj 2 \ 

M aP (uj 2 ) = a aP + b aP uj 2 + 0{lj 4 ) , (3.4) 

and hence Det \\M aP ( uj 2 ) — 5 a p\\ ~ A + Bu 2 . For the zero 
modes, A must vanish. In fact, it is not zero due to different 
approximations, such as, using a grid in (E, L) sub-space, 
substitution of the infinite matrix M aP by a matrix of finite 
size ( N a x N a ), and the replacement of an infinite series in 
l\ by finite series of length (fi) ma x. Assuming that an error 
is Slo = \A/B\ 1 ^ 2 , we can select the best basis by testing 
different functions g(r). It turns out that for Lt = 0.1 the 
appropriate choice is g(r) = —p(r), with 5uj ~ 0.03. The zero 
mode test allowed us to verify the accuracy of the equations 
and obtain the accuracy estimate of eigenfrequency calcu¬ 
lations. We also note that the accuracy drops sharply at 
Lt < 0.010 

For g(r) = p'o(r)/<&' 0 (r), the shear mode eigenfunction 
for the potential consists of only one element: x a ( r ) = 
df/r) J 1 ’ 11 . So, all the elements in the first row and first 
column of matrix M°‘ p [fi) must be close to zero, except that 
the first one is equal to unity. 

The results of the stability study are the following. 

1. There are no unstable solutions corresponding to odd 
spherical harmonics l. 

2. For even values of l we found only aperiodic unstable 
solutions, Reu> = 0, uj = i'y. 

3. The unstable models are found within the range Lt < 


2 For very small Lt, we have developed recently a special ap¬ 
proach which allows to investigate spectrum of eigenfrequencies 
even for almost pure radial models, Lt —> 0. This approach and 
its application will be presented in separate work. 
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Figure 3. Radial profiles of the potential perturbation for the 
model Lt = 0.01 corresponding to different modes of unstable 
spherical harmonic 1 = 2. To save room the profiles are shown 
in a single plot, but vertically separated from one another ( 7 ^ 
decreases from bottom to top). 

0.316, which corresponds to £ = 2T r /T± > 2.2, or global 
anisotropy £ > 0.55. 

4. For a given spherical harmonic l , the number of 
unstable modes increases infinitely with decreasing Lt . 
These modes have different growth rates 7 j 1 \Lt), j = 
1,2, ...,j$L(L T ) (see Fig.0. 

5. Eigenfunctions of the radial part of the potential per¬ 
turbation, corresponding to different modes j differ in the 
number of nodes (see Fig. Larger growth rates 7 , corre¬ 
spond to modes with fewer nodes. 

6 . With increasing Lt, the growth rates 7 ® decrease, 
and modes with large number of nodes disappear. 

Stabilization of all modes of l = 4 harmonic occurs at 
Lt = 0.142 (( = 5.1; £ = 0.8); all Z = 6 modes stabilize at 
Lt = 0.077 (( = 8.7; £ = 0.89). The largest-scale nodeless 
mode (Z = 2) remains the most unstable for any value of 
parameter Lt < 0.316. 

Fig. [2] compares the obtained growth rates of Z = 2,4, 6 
unstable modes, and some characteristic frequencies, such 
as radial frequency oscillations fir, and 


(fij)max = \/47rGp(0) , 



(Hj)loc = \J47 rGp(2ri) , 


where n is the isotropic radius of the nucleus. (flj) ave 
is a weighted Jeans frequency; (flj)LOC is a Jeans fre¬ 
quency on the radius of localization of the perturbation (i.e., 
tloc ~ 2 ri). It is seen that the growth rates are small com¬ 
pared with all characteristic frequencies (1331) : for instance, 
(Dj)loc exceeds the growth rates for more than an order of 
magnitude. 

The obtained slowness of the modes allows one to aver¬ 
age over the motion of a particle along its orbit, and consider 
orbital slow dynamics rather than ordinary particle dynam¬ 
ics (see, e.g., Polyachenko, 2004, 2005). Then the matrix 
equation can be simplified, and instead of the full equation 
based on matrix (13T2T) . one may consider a ‘slow’ equation, 
which is obtained from (13721) by omitting all terms except 


those for which Zi = —Z 2 / 2 . Similar simplification is used by 
Palmer (1994) for calculation and interpretation of the insta¬ 
bility in the ‘Lynden-Bell’ approach. Recall that for highly 
eccentric orbits Q .2 ~ fii/ 2 , which means 

ZlfZl + I2LI2 = I2 (fl 2 — | Ul) = I2 flpr <C Ui,2. ( 3 - 6 ) 

We checked the applicability of the ‘slow’ approach by direct 
recalculation of the spectra of modes for the lower spherical 
harmonics (Z ^ 6 ). The comparison demonstrates that dif¬ 
ference in frequency values does not usually exceeds 1 per 
cent. Thus we conclude that for q = 0 series ROI can be 
interpreted in terms of ‘Lynden-Bcll’ mechanism. 

This result seems suspicious in the absence of a dom¬ 
inant external potential, which provides a slow precession 
for all orbits and ensures the slow mode (Polyachenko et 
ah, 2010). An order-of-magnitude estimate for eigenfrequen- 
cies gives the Jeans frequency, or the dynamical frequency, 
oj ~ Qj ~ fldyn ~ \JGM/R 3 = 1. However, this is a hasty 
conclusion: although the characteristic dynamical frequency 
is of the order unity, the maximum dynamical frequency of 
star oscillations ‘locked’ near the centre is very high. For 
q = 0 series a large group of stars never leaves the cen¬ 
tre, and despite the obtained growth rates are substantially 
greater than unity, they are still much smaller than the dy¬ 
namical frequency of the locked stars. The slowness occurs 
here due to the small deviation of the potential from a har¬ 
monic form that exists in the central region. Note that these 
modes are turn out to be analogous to slow modes in near- 
Keplerian systems (Tremaine, 2001): in both cases orbital 
precession rates are low, and the modes are formed due to 
orbit-orbit alignment. 

We performed an ‘experiment’ to determine which par¬ 
ticles give the main contribution to the growth rates, retain¬ 
ing only contribution from E c to some E = _E max < 0 in the 
integrals over the horizontal line L = Lt, E c < E < 0. Our 
calculations confirm that the main contribution to the ma¬ 
trix elements comes from particles with apocentric distances 
t max(77, L = Lt) much smaller than unity. 

.Jeans instability mechanism suggests a simple depen¬ 
dence of critical parameter (Lt) crit, at which the system 
becomes stable, from the spherical harmonics Z. For a system 
in equilibrium oy ~ R y/Gp, where oy is the radial disper¬ 
sion. On the other hand, assuming marginal stability one can 
obtain for the Jeans characteristic scale in the transverse di¬ 
rection, A j, from the Jeans criterion A j ~ cr±/\/G p. Thus, 
using Aj ~ R/l, we have 

ih ~ IF ~ ’ (3 ' 7) 

which gives Lrit ~ (1 — £) _p with the exponent p of order 
unity. Since for small Lt, Lt — (1 —£)(l/2 — q) (Polyachenko 
et al., 2013), it is natural to expect (Lt) crit to be inversely 
proportional to some power of l: (Lt) crit oc Z _1,/p . 

Fig. 0 shows the stability boundaries (Lt)ci it for even 
harmonic numbers Z in the range 2 ^ l Sj 30. The filled cir¬ 
cles show the results obtained using the full matrix equation 
(HOI) . Starting from Z = 10, the linear combination of orbital 
frequencies = I 1 LI 1 +I 2 LI 2 vanishes for some Zi, Z 2 . Due 
to these resonances, calculation of the stability boundaries 
becomes extremely time-consuming. Open circles show the 
results of calculations using the ‘slow’ equation, which are 
almost identical in the absence of resonances (Z ^ 8 ). How- 
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Figure 4. The dependence of the critical parameter (Z/T)crit on 
the harmonic number l: q = 0 series (filled circles show the stabil¬ 
ity boundary inferred from the full integral equation EU , open 
circles from the ‘slow’ integral equation); g = — 1 series is shown 
by triangles (full equation). The dash-dotted line shows the least 
square linear fit to the q = — 1 boundaries. 


ever, when resonances appear, the results begin to diverge. 
The full equation gives approximately exponential decay for 
(Lt) ait (/), while the ‘slow’ solution decreases significantly 
faster than exponent. 

We conclude that the estimate for (Lr)crit(O based on 
usual Jeans relations for gravitating medium is incorrect to 
describe the bar-forming instability in highly heterogeneous 
systems. It enables only to predict the decrease with l. In 
addition, ‘slow’ solution is applicable only in the absence 
of resonances. Correct calculations of spherical harmonics 
l ^ 10 are possible by using the full equation only. 

Our matrix calculations was also supported by numeri¬ 
cal IV-body simulations using Superbox-10 code (Bien et al. 
2013). This code is an example of a particle - mesh scheme, 
which solves the Poisson equation by fast Fourier transform. 
The number of grid points N g for each coordinate is the same 
and is taken so that the number of grid cells be compara¬ 
ble with the number of particles. The number of particles 
in all calculations except one was 10 6 ; with N g = 256. The 
model close to the stability limit (Lt = 0.2975) was calcu¬ 
lated with 10 7 particles and N g = 512. The code uses three 
meshes. The biggest one allows us to simulate interaction 
between galaxies. Medium meshes are designed to simulate 
separate galaxies, and the smallest meshes are used to re¬ 
solve fine structures in galactic centres. In all our calcula¬ 
tions, the mesh sizes were taken to be 30, 5 and 1 (recall that 
the initial radius of the system R = 1). The growth rates, 
evaluated from numerical experiments show good agreement 
with results of matrix calculations, especially for the mod¬ 
els with moderate growth rates. On the other hand, in the 
models with Lt < 0.14 it is difficult to distinguish the tem¬ 
poral interval of exponential growth of perturbations. This 
can be explained by interference of l = 2 and l = 4 modes. 
To investigate these models, initial states of high degree of 
symmetry as well as filtering of higher harmonics of force 
are needed. 


3.2 q = —1 series 


In contrast to q = 0 series discussed above, models of q = — 1 
series are monoenergetic, e.g. all stars have the same energy. 
In the limit of purely radial orbits all of them reach the 
outer radius; thus stars locked near the centre are absent. 
Accordingly, the dynamical frequencies of stars are almost 
identical and are of the order unity. Here we investigate how 
this affects the stability properties. 

For the given q, the distribution function CD reduces 
to: 


F(E,L) 


N(—l, L t ) H(L t - L) 

47t 3 Lj. 


6(E), 


(3.8) 


where 6(x) is the Dirac J-function. In the limiting case of 
purely radial orbits, Lt = 0, we obtain a model discussed 
by Agekyan (1962). On the other hand, systems become 
isotropic at Lt ^ 0.5613. 

In the calculation of the matrix element, M aP (ui), two- 
dimensional integration is reduced to a one-dimensional in¬ 
tegration over the vertical interval E = 0, 0 < L ^ Lt- 
However, the integrand includes a derivative with respect 
to energy, so some functions should be found in a small 
neighbourhood A E < E < 0,0 < L < L c i IC (E). The final 
expression of the matrix element obtained from (tOl) has a 
form: 


M aP (uj) 


2 N _d_ 
Zf dE 


j LdLS a \E,L) 

_0 


E =0 


4./V 

Lt 



x 



d L f LdL 2 (E,L) 

dE J L t Qi(E, L)[uj 2 — Llf t (E, L)\ 
0 1 2 


(0, L T ) ipifiz (0, L t ) \ 
+ Qi(0, L T ) [u 2 - nf ih (0, L t )] j 


(3.9) 


where S aP (E,L) is defined by (13.31) . The first term on the 
r.h.s. of can be converted to 


2 N d 
If dE 


/ 


LdLS aP (E,L) 


2 N 

7T Lj, 


Lo 


E =0 


/ 

J'l 


X“(r) x P ( r ) P dr 

y/2V(r) - L 2 T Ir 2 


1 

/ 


X“(r) X /3 ( r ) P dr 


(3.10) 


where in the first integral on the r.h.s. of (13.101) tt ,2 are zeros 
of the radicand in the denominator. 

Stability of individual harmonics is investigated by em¬ 
ploying equation dj with the matrix elements (EH) for 
1 ^ l ^ 30 and Lt > 10~ 3 . Note that for the model with 
Lt = 0.01, parameter £ = 2 T r /T± = 172 and the global 
anisotropy £ = 0.994. 

The spectrum of unstable modes for models of this series 
differs signfficantly from the spectrum of q = 0 series models. 
This applies to both even and odd spherical harmonics l. 

1. For even l, there is only one aperiodic unstable mode, 

1) ( l t) = *To°( L t ) 


and several oscillating unstable modes 

Uj l) (L T ) = uj( 1) (L t ) + i^ l) (L T ), j = 1,..., j, 
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q=1 series 



Figure 5. Same as in Fig. [2] for q = —1 series. Stabilization of 
the most unstable mode / = 2 occurs at L cr ~ 0.47. The shaded 
area of Qi variation is degenerated into a thin almost horizontal 
line. 

in which the growth rates decrease with real part of the 
frequency increasing, 

< of(Lr),7j!!i(£T) > lf{L T ), 

and the number of unstable oscillatory modes j max depends 
on Lt- The real parts of the frequencies are separated by 
Hi (see Fig. [6] a): 

-uf± 1 nQ 1 (E = 0,L = L t ). 

2. For odd modes starting with l = 1, there are only 
oscillating unstable modes, the real part of the frequencies 
are approximately equally spaced (see Fig.[6]b). 

3. Growth rates of all modes decrease with increasing 
Lt- For a given even spherical harmonics l, less unstable os¬ 
cillatory modes stabilize first, then more unstable aperiodic 
modes stabilize. 

All non-spherical harmonics are fully stabilized when 
Lt > (Lt) crit ~ 0.47, corresponding to £ < 1.5 or £ < 
0.34. It coincides with the stabilization of aperiodic bar¬ 
mode instability, l = 2 (see Fig. 0. Stabilization of l = 4 
harmonic occurs at Lt = 0.284 (£ = 3.7; £ = 0.73) and 1 = 6 
harmonic - at Lt = 0.182 (£ = 6.9; £ = 0.85). 

Fig. [4] shows the stability boundary (Lt) crit for even 
modes in the range 2^1^ 30. The unstable modes (shown 
by triangles) are obtained using the full matrix equation 
(13791) . They fit well a simple relation: 

ln(L T )crit = —0.2359 1 — 0.2717 . 

Slight deviation from linearity at l ^ 26 may be due to 
insufficient accuracy of the calculations. 

Fig-El shows oscillatory and aperiodic modes of even 
harmonics for different values of Lt- Closed symbols show 
the solutions obtained using the full equation (BUl) , whereas 
open symbols indicate eigenmodes found with a simplified 
equation derived from (El) by neglecting the terms asso¬ 
ciated with the energy derivative of the DF. The remain¬ 
ing last term in am mostly determines the eigenfrequen- 
cies, which also follows from the figure. For strongly radially 
anisotropic models, the growth rates of aperiodic modes are 



0 1 2 3 4 5 

Re co/ll-i 


Figure 6. The oscillating modes of q = —1 series for three values 
of Lt ( Lt = 0.01; 0.1; 0.3): (a) l = 2, (b) l = 3. Solution (E) 
obtained by exact matrix equation (l3~9l) : solution (A) by an ap¬ 
proximate equation (see explanation in the text). All frequencies 
are given in units of the radial frequency for purely radial orbits 
Hi « 2.16. 

much larger than unity. However, at Lt > 0.1 they become 
of the order unity and comparable with the growth rates of 
oscillatory modes. Note that the radial frequency Hi (E,L) 
enters the equation at the boundary term E = 0, L ^ Lt 
only, and for small Lt is close to a limit Hi (0,0) « 2.16. 
Thus, all the modes in this series can not be considered as 
slow ones. 

Note that in contrast to q = 0 models, in which some 
stars are locked in the centre, in q = — 1 models the orbits 
of stars have large radial excursions and are no longer in the 
regime of slow dynamics. This is a possible reason for the 
fastness of the q = — 1 modes. 


4 CONCLUSION 

This paper analyses two approaches presented in the liter¬ 
ature to interpretation of the radial orbit instability (ROI). 
The first one explains ROI in terms of the classical Jeans 
instability. Indeed, it is natural to expect that radially 
anisotropic systems are ‘cold’ enough in the transverse di¬ 
rection for the instability to develop. Such an approach has 
been proposed in the first works devoted to ROI (Poly- 
achenko and Shukhman 1972, 1981). The second approach 
appeals to a bar formation in disc galaxies proposed by 
Lynden-Bcll (1979). His mechanism considers coalescence of 
the so-called ‘abnormal’ stellar orbits, for which the preces¬ 
sion rates of stars, H pr , decreases with decrease of angular 
momentum L while the adiabatic invariant Jf = Ir + | L is 
conserved (here Ir is the radial action), i.e. ( dLl pr /dL)j f > 
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0. Let a weak bar-like perturbation of the potential rotates 
at a slow rate Q p . Orbits having a form of nearly symmet¬ 
ric ovals and the precession rates close to bar pattern speed 
effectively interact with the bar. The torque exerted by the 
bar on the abnormal orbit results in a change of the orbit 
precession rate, so that the orbit tends to line up with the 
bar, contributing to the potential well, thus enhancing the 
bar-like perturbations. 

For spherical systems, highly eccentric orbits usually 
obey required inequality dfl pl /dL > 0. Therefore, the idea 
of Lynden-Bell can allegedly be extended to spherical sys¬ 
tems. However, the key assumption that makes this analogy 
legitimate is slowness of the perturbations, i.e. perturbation 
frequencies must be much smaller than the characteristic or¬ 
bital frequencies. Otherwise, the concept of orbit as a sepa¬ 
rate united object (instead of a set of individual stars) inter¬ 
acting with potential perturbation is invalid. Moreover, the 
orbital approach can be used to describe symmetric relative 
to the center disturbances only. The easiest way to explain 
it is to consider perturbations on a disc, which have a form 
5$ oc cos {rmp). Even m describe symmetric perturbations 
relative to the center, while for odd m signs of the potential 
are opposite at the opposite points: <f>(r, ip + 7r) = — >f>(r, i p). 
In the latter case, the potential exerts differently on either 
side of the symmetric elliptic orbit. Note that for spherical 
models, the role of the azimuthal number m plays a spher¬ 
ical number l. If an unstable model violates any of these 
assumptions (the perturbation is not slow, or instability is 
possible for odd l), then the validity of the ‘Lynden-Bell’ or, 
equivalently, ‘slow even-F or ‘orbital’) interpretation for this 
model can be put into question. 

In this paper, we consider two series of DFs of the form 
F(E,L). In the first one, 5 = 0, unstable perturbations are 
indeed slow, and are possible for even l only. In the second 
series, q = — 1 , the modes are fast, and instability is possi¬ 
ble both for even and odd l. The reason for this difference is 
in the properties of energy dependence of their DFs. In the 
model with equipartition of the energy (q = 0 ) many stars 
never leave the central region (their apocentric distances are 
much less than the radius of the system), where frequency 
of radial oscillations, fii, and Jeans frequency, flj, are much 
greater than frequencies of the perturbation. In this sense 
unstable modes are indeed slow, and thus the orbital ap¬ 
proach to the instability is valid. 

On the contrary, in the mono-energetic model (5 = — 1) 
orbits of all stars nearly approach the outer radius, so our 
estimates that mode frequencies would be comparable to 
the characteristic dynamical frequency (almost the same for 
all stars) is proved true. Hence there is no way to use the 
orbital approach, and the odd modes obeying |(u| > fli can 
be excited. 

We conclude that the spectra of radially-anisotropic 
DFs of the form F(E,L) are allowed to have both slow 
and fast modes. However, some DFs support essentially slow 
modes, while others allowing for both slow and fast modes. 
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APPENDIX A: LAGRANGIAN FORM OF 
MATRIX EQUATION FOR SPHERICAL 
SYSTEMS 

Matrix elements in equation (EH) normally has the form 


M a \u) = A-kG{2'k) a 


ii = — 0012 = — l 


dELcLL 

o 


n hh (E,L)Z | + z 2 ^ 

<^ E ’Q U - h *™E,L) dL ’ (M) 


which is sometimes called ‘Euler’ form of the matrix ele¬ 
ments. The matrix contains E- and L-derivatives of the DF 
F(E,L). Our notations are explained in the main text. 

Integration is taken over domain T> of the two- 
dimensional sub-space ( E,L ), which depends on the DF. 
Boundaries of T> consist of line of circular orbits £ C irc, line 
of radial orbits £ r . o., and a line of escape ^singular, which is 
usually given by E = 0. 

The Euler form of the matrix element has a significant 
drawback: if the DF F(E, L) is singular on any of the bound¬ 
ary lines (but the singularity is integrable), the Euler form 
contains divergent integrals. Note that two natural phase do¬ 
main boundaries £ c \ic and £ r . 0 . are safe in this sense, because 
there is no flux through these lines (see below). 

Our aim is to rewrite the Euler expression for the matrix 
element in the form that is free of E- and L-derivatives of 
the DF. For this we first write m as follows 


M aP {uj) 

= 47tG(27t) 2 



OO l C C 

£ / 


dELdLtff h (E,L) 


1 1 = —OO ln= — l 


UJ — 1 1^1 — ^ 2^2 




d ( l 


+ dL V fii F 


(A2) 


According to the Gauss theorem 

J I dEdLdiv^Ai^Bi^) = J d£(Ai 1 i 2 Bi 1 i 2 ) , 

e 

where £ denotes the boundary of T>, 

£ f-circ + £i.o + ^singular - 

The DF F(E, L) is allowed to be finite or to have integrable 
singularity on Singular- 

Vector A ■ B vanishes on the boundary £ T . 0 ., as well 
as the vector flux across the boundary. This is due to the 
factor L in the expression for Bi ± i 2 . Furthermore, flux across 
the boundary £ C ii C vanishes due to the identity (E, L ) = 
<$iiOX“(-R), i.e. it vanishes on the circular orbits for all h and 
1 2 , except li = 0. 

Yet, flux across Singular is not zero, if DF takes finite 
values at this boundary. Moreover, the flux is infinite if the 
DF is singular on Singular. To summarize, 


j di ( AB ) = 0 , J d£ ( AB ) = 0 , J d£ ( AB ) ^ 0 . 


^r.o. 


Thus, the Euler form of the matrix element is invalid 
when the DF is singular but integrable on one of the bound¬ 
ary lines. Besides, non-integrable singularity is present in 
the integral EH- We may conclude that standard lineariza¬ 
tion procedure is incorrect near the boundaries if one uses 
the Euler technique. It is therefore clear that the attempt 
to bring the Euler expression to valid (Lagrangian) form 
without derivatives of the DF by using the Gauss theorem 
(roughly, using integration by parts) was bound to fail. 

The desired Lagrangian form could be obtained if the 
flux through £ vanished. In this case, we have 

Jhh = ~J J dEdL (Ai 1 i 2 S7)Bi 1 i 2 , 


using the identity 


JL0l +JL_L = 0 

dE Lh + dLQ 1 


Now we introduce a ‘vector’ A tl i 2 = [[Ae)i 1 i 2 ,{A l )i 1 i 2 S J 
with coordinates 


0,2 


and denote 


Bi li2 (E,L) = 


h 


(Ae)^^ = ( Zi + I 2 ^ ) F , ( A L ) h i 2 = — F 


fil 


L^{E,L) 

uj — liLli — I 2 LI 2 


Then the integral 

dELdLtff h {E,L) 


Jhh= .[ J 


_d_ 

dE 


UJ — li^li — I 2 LI 2 

( , ‘ +fc t) F 


+ A(^f 


dL V fii 


(A3) 


can be written in the form: 


dEdLBi x i 2 divAqi 2 = 


Jhh ~ 11 

J J dEdL ^div(Ai li2 Rq !2 ) - (Aip 2 V)Rq; 2 j . (A4) 


and for the matrix element: 


1 


dEdL 




M aP (w) = -AnG(2n) 2 £ £ D\ 

l\ = — 00 l 2 = — l 

x i Qhh dE +12 


f(e,l ) 


L <u 


■)• (A5) 


However, to obtain this (correct) form, one needs to use 
the Lagrangian technique from the very beginning. It was 
first presented in the works by Kalnajs (see, e.g. Kalnajs, 
1977). Rewriting a matrix element for disks (equation (16) 
of his paper) to the spherical geometry (in which we are 
interested in) one can have for the perturbations which are 
independent on the azimuthal angle 


M ap {w) = -47rG 


21 + 1 


??/// dhdl2dI 3 F{E,L) 


f, d d 

V 1 dh +l2 dh 


A>rf, 2 (h,i2,h) 


w — n lll 2 (/i,/ a ) 


(A6) 


Here Ii is the radial action, J 2 = L, I 3 = L z , fip 2 = 
dE{I 1 ,I 2 )/dI 1 , 2 , = 0 . Expression h (I u I 2 ,1 3 j de¬ 

notes 


Kiz ( 7 i- ^ 7 3) = i2 (/ 1 , h, h) i2 {h,h, h )]*, (A7) 
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where <!?/[ i 2 (Ii, I 2 , 13 ) can be written in the form (see also 
Appendices in Polyachenko et al., 2007, Polyachenko and 
Shukhman, 1981) 

$? li2 (E, L,0 O ) = 2tt P / 2 (0)Pf * 2 (sin 0 O ) ( E , L) , 

(A 8 ) 

where sin#o = L z /L, — | n < 60 < 1 n. From (1 ATI) and (IA 8 I) 
one finds 

= (2-K) 2 P^(0)P l - h (smdo)e ih ^T 1 i 2 (E,L) 

x P- 1 - (0)P / 2 (sin O 0 ) e-^4,\ t _ l2 (E, L ) . (A9) 

or 


= (2-7t) 2 P[ 2 (0) P i _, 2 (sin 6 | o) 


x Pi~ l2 (0) P} 2 (sin # 0 ) ip“fi (E, L) . (A10) 


hh 

Then one needs to change the differentiation with respect to 
actions 7i and I 2 to differentiation with respect to E, L, do 
in (1A6I) : 

d 


, d , d d , 

h -dh +h dh =^dE +h 


d sin d 0 


dL L i9(sin6>o)_ 
transform the volume element dI\dI 2 dIo = 
dE L dL d(sin do )/fh, and integrate over z = sinPo. 
Since P/ 2 (± 1) = 0 provided I 2 7 ^ 0, and 

£ dz z d[P / 2 (z) Pr h (z)] dz = -J Q 1 dzP}* (z) Pr h (z), 
then 


hi 


d z d 
~dL T~dz 


KU E,L,z) 


1 

= {jL + l) f dz ^ {E ’ L ’ z) ■ (A11) 

-1 

itegral over z is 
£ dzPt*(0)Pr l2 (z)Pr l2 {0)Pl 2 (z) = [2/(2 1 + 1)] D[\ 


Given that the integral over z is 


i-i 

one finds 
1 


/ Kl (E,L, z) dz = (27r) 2 D \* {E ,L ) . (A 12 ) 

-1 

Substituting (1A12I) into the r.h.s. of (lAllf) . and then substi¬ 
tuting the resulting expression in (IA6I) . we obtain 




dELdL 


F(E,L) 


dE 


+ h 


JL I 

dL + L 


Hi 

fdzVtf 2 (E,L,z) 


IV — 


(A13) 

Finally, we have the desired expression of the matrix element 
in the Lagrangian form: 


OO l 


G(2„) 2 £ T. D ‘‘ f 

7 . __ 7 _ 7 J J 1 


1 1 = —OO l'2 = — l 

d . , d\fL1>?X 


■ (A14) 


APPENDIX B: CONSTRUCTION OF THE 
BIORTONORMAL BASIS SETS 


The effectiveness of the matrix method m depends on 
the proper choice of a basis function set {x a ( r ): P a ( r )}> a = 
1 , 2 ,... satisfying the requirement of orthogonality: 

(x“ P P ) = f x“0) p P {r)r 2 dr = -5 ap (Bl) 

JO 

and the Poisson equation: 
d 2 


2 d 1(1 + 1) 
dr 2 r dr r 2 


X a (r) = p“(r) . (B2) 


This basis is used in expansion of the perturbation potential 
and density (1231) . 

Instead of (IB2I) . let us consider the following eigenvalue 
problem, assuming l ^ 1 : 

X(r) = Xg(r)x(r) , (B3) 


1(1 + 1) 


_ IjL 

dr 2 r dr r 2 
with boundary conditions 

x(0) = 0, x'(l) + (l + l)x(l) = 0 ■ (B4) 


The solution is a discrete set of positive eigenvalues, A“, and 
eigenfunctions, X a ( r )i orthogonal with weight r 2 g(r): 

1 

J drr 2 g(r)x a (r)x^( r ) 0C d a ^ . (B5) 

0 

Then the biortonormal set consists of functions x“ (r) for the 
potential and functions p a = A“ g(r) X a ( r ) for the density. 
According to (IbsTi . the biortonormal condition m will be 
satisfied if the normalization of the functions x a ( r ) obeys 



(B6) 


In particular, for g(r ) = —1 we obtain the well-known 
biorthogonal basis: 


Xi ( r ) = 

P?(r) 


y/2 1 

x a |Ji + i/ 2 (>r a )| 

I-^+1/2 (^a) I 


Jl + l/2(Xa r) 

sh 

Jl+l/2(Xa r) 


(B7) 

(B 8 ) 


where yc a are the positive roots of J;_ i/ 2 (xa) = 0 , a = 
1,2,... (Polyachenko, Shukhman 1981). 

Choice of g(r) = p'o( r )/&o( r ) is interesting for calculat¬ 
ing the shear zero mode 1 = 1, with the perturbed potential 
and density, respectively, x( r ) = $o( r ), an d n(r) = p' 0 (r). 
The expansion for (12.51) contains only one term correspond¬ 
ing to the minimum eigenvalue A = A “ =1 = 4-77 (i.e., C 1 = 1; 
C k = 0 for k > 1) and the basis function x 1 ( r ) = 4?o( r )- 
Note that Bertin et al. (1994) discuss a similar technique 
for infinite systems. However, there is no attempt to build a 
system so that the expansion of the potential for shear mode 
contains a single function (which is useful for tests). 

The eigenvalue problem (IB3I) in the form of a differential 
equation can be reduced to a more convenient eigenvalue 
problem in the form of an integral equation. In doing so, we 
write the equation (lB3l) in the equivalent form: 

1 

X(r) = ~ 2l + l / dr ' r ' 2 9 ( r ')x(r')Fi(r,r') , (B9) 

0 





























On the radial orbit instability 


where J-(r, r') = r< = min(r, r'), r> = max(r, r'). 

In this form, the boundary conditions are satisfied automat¬ 
ically. 

Note that construction of a basis set for radial pertur¬ 
bations (l = 0) is not covered by above analysis, and should 
be discussed separately. 


